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SECTION  1 
INTRODUCTION 

It  is  well  known  in  aerodynamics  that  shockwave-boundary  layer  interac¬ 
tions  are  dependent  upon  local  viscous  dissipation  and  thermal  conduction 
processes  for  the  distribution  of  internal  and  kinetic  energies.  This  depen¬ 
dence  can  in  many  applications  be  a  sensitive  determinate  of  macroscopic  flow 
properties.  In  airblast  applications  it  is  clear  that,  in  both  clear  and 
dusty  atmospheres,  irregular  shock  reflection  processes  evolve  from  such 
local  shock  interactions  with  boundary  layers  and  possibly  with  other  inter¬ 
nal  shear  layer  regimes.  These  types  of  shock-boundary/shear  layer  processes 
impose  severe  demands  upon  computer  codes  which  attempt  to  calculate  airblast 
effects. 

It  is  certainly  important  in  studies  of  shock-boundary  layer  inter¬ 
actions  to  either  eliminate  or  suppress  to  extremely  low  levels  the  purely 
numerical  dissipation  effects  which  are  present  in  existing  hydrodynamics 
codes.  This  is  both  difficult  and  ambiguous  in  inviscid  PDE's  because: 
(i)  they  contain  no  physical  dissipation  operators  which  are  present  in  real- 
world  shocks  and  boundary  layers,  and  (ii)  the  codes  which  solve  such  PDE’s 
thus  generate  only  numerical  dissipation  in  those  critical  interaction 
regimes  of  the  PDE  solutions.  It  follows  that  the  most  reliable  computations 
of  the  essential  physics  of  shock-boundary  layer  interaction  processes  must 
first  include  the  physical  dissipation  operators  of  the  full  Navier-Stokes 
equations  and  then  use  optimally  distributed  grid  densities  which  can 
accurately  resolve  the  actual  physical  dissipation  processes  and  eliminate 
the  nunerical  dissipation  processes. 

It  appears  that  neither  of  these  last-mentioned  factors  have  been 
accommodated  adequately  in  existing  airblast  codes;  and,  for  fundamental 
physical  and  mathematical  reasons  there  is  cause  for  concern  that  these 
factors  may  never  be  handled  adequately  by  extended  versions  of  most  existing 
airblast  codes. 


As  a  possible  alternative,  the  Moving  Finite  Element  (MFE)  method,  which 
was  discovered  by  Professor  Keith  Miller  in  the  mid-1970's  has  shown  some 
promise  of  adequately  dealing  with  the  combined  physical  and  numerical 


requirements  of  irregular  shock  reflection  processes.  It  has  so  far  been 
shown  in  1-D  calculations  that  the  MFE  method  successfully  resolves  shock- 
boundary/shear  layer  interactions  according  to  the  full  Navier-Stokes  equa¬ 
tions.  In  these  1-D  results,  the  essential  physical  dissipation  processes 
were  resolved  accurately  over  highly  disparate  spatial  and  temporal  scales 
with  a  minimal  number  of  optimally  distributed  moving  grid  nodes.  But  despite 
the  promise  shown  in  these  early  1-D  results,  it  was  by  no  means  certain  that 
the  MFE  method  would  also  extend  effectively  to  2-D  for  the  resolution  of 
irregular  shock  reflection  processes. 

In  view  of  these  circumstances,  the  present  research  and  development 
effort  was  undertaken  by  DNA  as  a  two-year  project  to  determine  the  feasibil¬ 
ity  of  applying  the  MFE  method  to  airblast  problems  in  2-D  and,  if  warranted, 
to  develop  a  2-D  MFE  production  code  for  airblast  applications.  This  report 
summarizes  the  results  of  the  first  full  year's  effort  which  has  been  devoted 
to  the  following  work  statement  tasks: 

1.  Survey  and  select  sample  problems  which  test  the  effectiveness  of 
the  MFE  method  and  alternative  PDE  solution  methods  for  resolving 
blast  wave  effects,  including  shock  reflections  and  interactions  in 
2-D,  with  particular  attention  to  the  resolution  of  physical  vs. 
numerical  dissipation  effects. 

2.  Develop  new  coding  of  the  gas  dynamics  operators  in  2-D  for  MFE 
applications  to  hydrodynamics. 

3.  Test  alternative  node  control  strategies  for  2-D  hydrodynamics 
appl ications. 

4.  Test  and  implement  more  efficient  matrix  solution  methods  in  the 
Interest  of  computational  economy  in  the  MFE  hydrodynamics  code. 

These  tasks  have  been  carried  out  successfully  as  will  be  described  in 
the  three  main  report  sections  which  follow  sequentially  under  the  headings 
of  MAIN  TEXT,  CONCLUSIONS,  and  RECOMMENDATIONS. 


SECTION  2 
BACKGROUND 


The  problem  of  primary  concern  Is  the  accurate  calculation  of  HOB  over¬ 
pressures  as  a  function  of  range,  altitude,  and  time  for  various  field  condi¬ 
tions.  Despite  nunerous  ongoing  efforts  to  calculate  such  HOB  data  by  alter¬ 
native  computational  approaches,  significant  discrepancies  remain  between  the 
various  calculations  in  comparisons  against  each  other  and  in  comparisons 
against  experimental  data.  While  debate  and  analysis  of  these  respective 
discrepancies  continue,  there  remains  the  gnawing  possibility  that  certain 
essential  physics  processes  In  regions  of  large  gradients  which  drive  the 
real-world  systems  are  neither  Included  nor  otherwise  resolved  adequately  in 
existing  code  calculations.  It  Is  also  possible  that  some  of  the  essential 
driving  physical  processes  are  not  measured  directly  or  otherwise  distin¬ 
guished  by  detailed  cause-and-ef feet  experiments  in  field  measurements. 

This  project  thus  attempts  to  obtain  accurate  MFE  solutions  of  the  full 
Navler-Stokes  equations  for  Irregular  shock  reflection  processes.  It  appears 
that  the  highly  detailed  reflected  shock  data  obtained  in  the  laboratory  by 
Glass  and  co-workers1  provides  the  most  appropriate  data  base  for  bench¬ 
marking  codes  which  would  attempt  to  calculate  shock-boundary/ shear  layer 
Interactions.  We  concur  with  the  (paraphrased)  view  expressed  by  Glass  that 
codes  which  are  to  be  used  to  calculate  HOB  effects  in  non-ideal  atmospheres 
should  certainly  be  able  to  solve  accurately,  as  a  first  pre-requisite,  those 
physical  processes  which  account  for  the  observed  isopycnlc  data  in  carefully 
controlled  laboratory  experiments  of  regular  and  irregular  shock  reflections 
In  Idealized  atmospheres.  We  therefore  consider  the  calculation  of  Glass's 
shock  reflection  data  as  an  essential  set  of  standard  shock  test  examples  for 
DNA  airblast  problems. 

2.1  MATHEMATICAL  BACKGROUND  OF  THE  MFE  METHOD. 

In  order  to  present  a  coherent  report  of  work  to  the  broadest  range  of 
readers,  we  present  here  a  brief  mathematical  sketch  of  the  MFE  method.  Con¬ 
sider  a  general  system  of  partial  differential  equations  (PDE's),  U  =  L(U)  or 


•  «4*‘J 


AM 


Using  piecewise  linear  approxminants  of  u;  ...  uN,*  which  are  of  the 
form  u  s  u*  ♦  ny  ♦  p,  on  a  hexagonal ly  connected  triangular  mesh  (Figures 
1-3),  application  of  the  chain  rule  to  the  differentiation  of  uk  gives: 


u.  *£ak.ak^  +  x.gk^  ♦  y.yk^  .  where 
j  J  J  J 


odt^  * 


.  &kJ  ,  3Uk  .  Ykj  .  3uk 

^7  ’  ^7  ’  ^  ■ 


The  key  aspect  of  the  MFE  method  lies  In  the  rigorous  retention  of  the 
nodal  motion  terms  Involving  x  and  y.  In  most  conventional  PDE  solution 
methods,  the  nodal  motion  terms  are  either  set  to  zero  (fixed  nodes)  or 
allowed  to  assume  some  arbitrary  selected  values  (e.g.,  mean  fluid  velocities 
In  Lagranglan  codes). 

Alternatively,  the  MFE  method  uses  a  very  fundamental  criterion  for  the 

evaluation  of  grid  node  motions  (x  and  y)  at  every  time  step.  That  Is,  the 

basic  MFE  method  Is  formulated  by  requfrfng  that  all  of  the  time  derivatives 
•  •  •  • 

alj .  aNj,  Xj,  yj  be  found  at  each  instant  In  such  a  way  as  to  minimize 

the  norm  of  the  PDE  residual,  ||u  -  L { U )  ||  .  This  has  been  found  to 
minimize  PDE  solution  errors  and  to  dramatically  reduce  the  number  of  grid 
nodes.  In  order  to  handle  the  degeneracies  which  can  sometimes  occur  In  such 
a  parametrlzation  of  PDE  solutions,  regul arizatlon  techniques  are  applied. 
That  Is,  the  more  general  function, 

*  =  wj  ||  Uj  -  L1(U)  ||  2  ♦  ...  w2  ||  uN  -  Ln(U)  ||  2 

+  trlSgle  (ejhj  -  Sj)2  •  (3) 

al titudes 


*The  present  mesh  triangulation  has  been  chosen  for  simplicity.  Many  other 
grid  meshing  schemes  are,  of  course,  possible;  and  many  of  the  alternative 
schemes  should  be  Investigated  in  future  work.  It  Is  interesting  to  note, 
however,  that  nunerous  MFE  results  confirm  that  PDE  solution  approximants  of 
higher  degree  are  not  nearly  so  critical  for  attaining  high  levels  of 
stability,  accuracy,  and  resolution  In  an  optimal  moving  node  solution  method 
as  they  are  In  fixed  node  or  less  optimal  adaptive  PDE  solution  methods. 
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calculations,  and  cylindrical  co-ordinates  are  to  be  used  for  HOB  calcula¬ 
tions.  Numerous  node  control  strategies  have  been  tested  with  alternative 
penalty  function  formulations  in  the  regularization  terms  of  Equations(4) . 
The  treatment  of  boundary  nodes  Is  particularly  significant  and  suggests  that 
alternative  grid  mesh  trlangul arizations  are  needed  for  the  most  effective 
solution  of  Irregular  shock  reflection  problems.  Finally,  alternative  matrix 
solvers  have  been  tested  and  Implemented  In  the  stiff  ODE  solver  of  the  MFE 
alrblast  code.  Widely  varying  convergence  properties  and  computational  effi¬ 
ciencies  were  observed  for  various  candidate  matrix  solution  methods,  and  it 
appears  that  either  ADI  methods  or  our  recently  developed  iterative  multi  grid 
type  solver  will  be  most  desirable  for  Irregular  shock  reflection  problems. 


2.2  DISCUSSION  OF  RESULTS. 

The  most  significant  results  which  emerged  from  the  above  tasks  are 
highlighted  here.  (A  detailed  description  of  every  Important  code  Implemen¬ 
tation  and  every  important  code  test  performed  In  this  work  would  be  too 
voluminous  to  be  very  useful.) 


2.3  MOTIVATION  FOR  MFE  SOLUTIONS  OF  NAVIER- STOKES  EQUATIONS  FOR  SHOCKS. 

It  Is  useful  to  first  consider  the  development  of  the  Navier-Stokes 
equations  In  i-D  In  the  context  of  kinetic  theory.  Recall  that  the  fluid 
equations  can  be  written  as: 

£.£<p,).0  (5) 


3(pv)  .  3  ,2.  .  3  j  <01  .  m  .  (2)  .  I 

~3t  3x  (p*  >  '  '  57  i"  *  »  *  •••) 


(6) 

(7) 


A  zero-th  approximation  of  the  kinetic  theory  for  gases  uses  the  consti¬ 
tutive  relations  for  an  Invlscld,  non-conducting  fluid;  l.e.. 
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p  =  p 


=  -tX-=  (y  -  1HE  -^pvfc) 


(8) 


q  =  q{0)  =  0  (9) 

These  relationships  yield  the  well-known,  inviscid  Euler  equations. 


A  first  approximation  of  the  kinetic  theory  for  gases  gives  the  Navier- 
Stokes  equations  according  to: 


p  =  P(0)  +  P(1)  =  (y  -  1) (E  -  Jpv2)  +  |y|£ 


(10) 


q  = 


(ID 


A  second  approximation  of  kinetic  theory  gives  the  Burnett  equations, 

etc. 
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The  issue  of  primary  concern  centers  around  the  potential  significance 
of  the  physical  dissipation  processes  as  determinates  of  macroscopic  flow 
properties  In  transient  shock  interaction  regimes.  From  the  standpoint  of 
physics  computations,  we  must  distinguish  the  relative  roles  in  code  calcu¬ 
lations  of  numerical  dissipation  vis  a  vis  the  physical  dissipation  effects 
which  appear  on  the  right  hand  sides  of  Equations  (6)  and  (7).  A  simple 
shock-boundary  layer  reflection  process  has  recently  been  solved  by  the  MFE 
method  in  1-D,  which  addresses  this  central  issue  directly;  it  is  the  so- 
called  anomalous  wall  heating  problem.  In  this  problem,  anomalously  high 
temperatures  are  calculated  by  existing  shock  hydrodynamics  codes  for  the 
reflection  of  a  planar  shock  in  an  ideal  gas  from  an  infinitely  reflecting 
wall  in  slab  geometry  (or  for  an  imploding  shock  reflecting  from  the  origin 
in  spherical  or  cylindrical  co-ordinates).  The  MFE  results  which  follow 
Indicate  that  the  anomalous  aspects  are  eliminated  when  numerical  dissipation 
Is  suppressed  and,  most  Important,  when  the  physical  dissipation  processes  in 
the  Navier-Stokes  equations  are  accurately  resolved  in  the  transient  reflec¬ 
tion  process.  The  following  test  problem  in  1-D  slab  geometry  illustrates 
these  results: 
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Initial  Conditions: 

p(  x ,  0)  =  1 
p(x,  0)  =  e(x,  0)  *  0 
v(x,  0)  =  -1 
v(x,  0)  =  linear 
v(x,  0)  *  0 
Y  =  5/3 


0.  <  x  <  2. 
0.  _<  x  <_  2.0 
Ax0  <  x  <  2.0 
0  _<  x  <  AX0 
x  s  0. 


Boundary  Conditions: 

Reflection  at  x  *  0. 

Dirichlet  at  Initial  values  at  x  5  2.0 


Ranklne-Hugoniot  Solutions 

s  =  1/3 

P+  =  4.0  ; 

e+  =  0.5  ; 

v+  =  0.  ; 

p+  «  1.33  ; 


for  Infinite  Shock  (t-*-00): 


p-  ■  1.0 

e-  =  0. 
v-  =  -1. 
p"  3  0. 


The  time  evolution  of  this  shock  was  solved  by  the  MFE  method  In  two 
ways:  First  the  full  Navler-Stokes  equations  were  solved  accurately  using 
alternative  values  of  v  *  4y/3  and  k.  These  solutions  are  denoted  by  N-S  In 
the  accompanying  figures.  In  one  set  of  N-S  solutions,  v  *  k  3  0.01,  which 
Is  unrealistically  large  but  which  permits  comparisons  to  other  fixed  node 
PDE  solutions  that  may  use  on  the  order  of  100  to  several  hundred  grid  nodes. 
In  another  set  of  N-S  solutions,  v  =  k  -  0.001,  which  approximates  physically 
realistic  values  for  actual  dissipation  processes  in  gases.  Second,  the 
variable  e  denotes  Internal  energy  per  unit  mass  in  the  accompanying  figures, 
and  anomalous  diffusion  (denoted  by  A  in  the  accompanying  figures)  was  simu¬ 
lated  by  Including  a  diffusion  term  vr  pxx  in  Equation  (1).  This  effectively 
simulates  some  form  of  uncontrolled  numerical  diffusion  which  is  present 
Intrinsically  In  all  numerical  PDE  solution  methods,  and  is  the  only  source 
of  dissipation  In  solutions  of  the  Invlscld  Euler  equations.  In  the  results 
which  follow,  we  have  verified  that  the  MFE  solutions  of  the  Navler-Stokes 
equations  have  reduced  all  numerical  or  other  diffusion  effects  to  Impercep¬ 
tible  low  levels  and  that  the  observed  behavior  of  shock  interactions  is 
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associated  with  physical  dissipation  processes. 

At  t  *  0+,  the  shock  Incident  on  the  origin  is  In  the  Incipient  state  of 
outward  reflection.  At  t  *  0.05,  Figures  4  and  5  show  that  the  calculations 
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Figure  4.  MFE  solutions  of  the  Nav i er-Stokes  equations  and 
anomolously  dissipative  equations  at  time  -  0.050 


of  the  reflected  shock  with  uncontrolled  diffusion  (or  simulated  nunerical 
diffusion)  tend  immediately  to  overheat  in  e  and  to  correspondingly  under¬ 
shoot  in  p  relative  to  the  Navier-Stokes  solutions.  Although  these  transient 
solutions  are  not  near  their  steady  state  values  at  this  early  time,  it  will 
be  seen  that  the  ensuing  evolution  toward  equilibrium  is  quite  sensitive  to 
both  the  magnitudes  and  the  nature  of  the  dissipation  processes  in  the  compu¬ 
tations. 

Figure  6  shows  that,  at  t  *  0.15,  the  lip  of  the  shock  in  the  Navier- 
Stokes  solution  is  approaching  the  steady  state  value  of  p  =  4.0,  and  tne 
anomolous  dissipation  solution  lags  by  a  significant  margin.  The  fluid  build¬ 
up  at  the  front  of  the  shock  is  evident  here  because  the  fluid  near  the  ori¬ 
gin  has  stagnated  while  additional  fluid  continues  to  stream  in  toward  the 
origin  from  the  region  to  the  right  of  the  shock.  Figure  7  shows  that  the 
anomolous  dissipation  results  continue  to  lag  behind  the  Navier-Stokes  solu¬ 
tions  to  a  significant  degree  at  t  *  0.300.  At  t  =  2.0,  the  Navier-Stokes 
solutions  have  approached  steady  state  Rankine-Hugoniot  conditions  (not  shown 
in  Figure  8),  and  the  anomalous  dissipation  solution  has  still  not  reached 
the  Rankine-Hugoniot  values  in  the  vicinity  of  the  origin.  The  anomalous 
wall  heating  effects  due  to  uncontrolled  dissipation  in  the  density  equation 
have  thus  persisted  to  very  long  times  vis  a  vis  the  accurate  solutions  of 
the  Navier-Stokes  equations.  Non-physical  dissipation  which  is  sometimes 
Introduced  as  artificial  viscosity  terms  In  fluid  calculations  can  be  shown 
to  have  similar  anomalous  effects,  as  can  the  numerical  dissipation  and/or 
nimerical  uncertainties  which  are  present  in  solvers  of  Inviscld  Euler  equa¬ 
tions. 

Figures  9  and  10  present  the  results  of  another  test  of  sensitivity  of 
the  Navier-Stokes  solutions  to  non-optlmal  grid  locations.  In  this  test  pro¬ 
blem,  a  physically  realistic  value  of  v  =  0.0001  Is  used  in  MFE  solutions  of 
the  Navier-Stokes  equations.  We  have,  however,  deliberately  constrained  the 
MFE  grid  nodes  In  this  test  case  so  that  they  do  m>t  migrate  to  their  truly 
optimal  locations,  as  In  the  results  considered  previously.  Figure  9  shows 
several  significant  features:  (i)  the  shock  gradients  associated  with  v  = 
0.0001  are  extremely  large;  the  accurate  resolution  of  these  gradients  would 
require  several  thousand  nodes  if  a  fixed  node  POE  solution  method  were  to 
be  used,  (11)  the  Rankine-Hugoniot  solutions  are  approached  much  more  rapidly 


Figure  6.  MFE  solutions  of  the  Navier-Stokes  equations  and 
anomously  dissipative  equations  at  time  =  0.150. 


MFE  solutions  of  the  Navier-Stokes  equations  and 
anomously  dissipative  equations  at  time  =  0.300. 


Figure  8.  MFE  solutions  of  the  Navi er-Stokes  Equations  and 
anomously  dissipative  equations  at  time  =  2.0. 
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for  the  physically  realistic  values  of  v  than  for  the  larger  values  of  v  which 
are  typically  used  either  tacitly  or  explicitly  in  many  other  POE  solution 
methods,  and  (iii)  the  slight  constraint  on  node  movements  and  thus  on  nodal 
positions  do  not  show  up  immediately;  but  once  the  perturbation  becomes  signi¬ 
ficant  (as  seen  in  Figure  10),  its  effects  can  grow  rapidly.  In  summary,  these 
results  demonstrate  that  reflected  shock  solutions  can  be  very  sensitive  to 
non-physical  dissipation  effects  and  to  slight  deviations  from  optimal  grid 
node  positioning,  even  in  adaptive  gridding  methods.  All  of  the  results  in 
this  section  were  obtained  with  approximately  30  MFE  nodes.  As  many  as  61 
MFE  nodes  were  used  to  verify  that  the  MFE  solutions  were  in  fact  converged 
solutions.  We  have  also  compared  the  steady  state  MFE  solutions  of  fluid 
velocities  with  the  analytic  Navier-Stokes  solutions  of  the  steady  shock. 
Excellent  agreement  was  obtained  between  the  MFE  solution  and  the  analytic 
Navier-Stoke  solution.  In  contrast,  it  was  found  after  extensive  attempts 
that  computed  solutions  which  contained  nominal  degrees  of  numerical  or  other 
anomolous  diffusion  could  not  be  forced  by  parameter  manipulations  to  agree 
with  the  analytical  shock  solutions.  In  still  other  MFE  calculations,  we  have 
computed  in  1-D  the  interaction  of  a  shock  with  an  internal  shear  layer  asso¬ 
ciated  with  a  contact  discontinuity  —  in  direct  analogy  to  the  irregular 
shock  reflection  mechanisms  which  occur  in  higher  dimensions.  Here  also,  the 
computed  macroscopic  flow  properties  were  sensitive  to  accurate  resolution  of 
the  physical  dissipation  processes  in  the  full  Navier-Stokes  equations.  In 
the  absence  of  the  stringent  tests  of  convergence  which  were  applied  here, 
it  can  be  extremely  difficult  to  discern  physical  oscillations  and  dissipation 
effects  from  non-physical  and/or  purely  nunerical  oscillations  and  dissipation 
effects.  We  have,  in  fact,  found  it  generally  impossible  to  simulate  the 
effects  of  the  actual  physical  dissipation  processes  in  shock-boundary  layer 
interactions  by  the  use  of  either  artificial  or  parametrically  controlled 
numerical  diffusion  processes  in  Euler  equation  solutions. 

2.4  MFE  DEVELOPMENTS  IN  2-D. 

The  2-D  MFE  airblast  code  is  being  benchmarked  extensively  against  what 
we  consider  standard  test  problems.  A  first  benchmark  is  the  wall  heating 
problem  discussed  above,  and  a  second  set  of  benchmarks  is  the  shock-on-wedge 
data  of  Glass  and  co-workers  for  the  various  types  of  regular  and  irregular 
reflections  of  plane  shocks  in  Ar.  All  of  these  problems  are  addressed  in 
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Figure  11.  Mf£  initial  grid  zone  configurations. 


rectangular  coordinates,  as  shown  by  the  schematic  view  of  an  initial  grid 
mesh  configuration  in  Figure  11. 

It  was  found  in  this  work  that  the  mass,  momentum,  total  energy  repre¬ 
sentation  (p ,  m,  E)  of  the  fluid  equations  has  better  numerical  properties 
for  MFE  code  calculations  than  the  velocity  representation.  The  Navier-Stokes 
equations  in  the  (p,  m,  E)  representation  are  given  by 

~  +  v  '  m  =  0  (12) 


—  +  V  •  (m  m/p)  =  V  •  4  (13) 

“  ♦  V  *  (Em/p)  =  - V  *  Q  +  V  ‘  (5*  (m/P))  ,  (14) 


where  Zis  the  generalized  stress  tensor  and  ()  is  the  heat  flux  vector.  The 
heat  flux  vector  is  frequently  given  by 


where  T  is  the  temperature.  The  general  stress  tensor  is  frequently  given  by 


“Pj.  ♦  _T 


(16) 


where  J  is  the  identity  matrix.  For  an  ideal  gas  (a  Newtonian  fluid)  in 
Cartesian  co-ordinates 
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It  is  important  to  note  in  Equation  (13)  that  the  quantity  m  m  is  a 
dyadic.  Special  attention  has  been  devoted  to  exacting  evaluations  of  those 
inner  product  terms  which  involve  factors  of  —  for  small  p.  These  factors 

p 

are  present  in  all  standard  representations  of  the  fluid  equations.  The 
inner  products  In  the  equations  above  have  been  coded  for  calculations  in 
both  Cartesian  and  cylindrical  co-ordinates. 


The  inner  product  formulation  of  discretized  finite  element  equations 
offers  a  potentially  unique  advance  for  the  treatment  at  the  origin,  as  r  . 
0,  of  those  POE  operators  which  are  proportional  to  1/r.  Singularities  do 
not  occur  because  the  Inner  products  of  the  MFE  basis  functions  are  weighted 
by  the  measures  2mr  in  cylindrical  co-ordinates.  The  MFE  method  is  thus 
1 ntri nsical ly  analytic  at  the  origin  in  both  cylindrical  and  spherical 
coordinates  which  obviates  any  of  the  numerical  tricks  which  are  otherwise 
used  in  discretized  PDE  solution  methods.  It  has  been  tested  and  verified 
that  there  is,  of  course,  a  price  to  be  paid  for  the  exploitation  of  this 
analytic  feature.  It  lies  in  the  fact  that  node  movement  properties  differ 
from  one  to  another  weighting  scheme  jf  MFt  norms.  Because  the  MFl  method 


is  still  a  new  POE  solution  method  In  many  respects,  the  specific  node  move¬ 
ment  strategies  for  cylindrical  co-ordinates  will  have  to  be  tested  much 
further  when  HOB  calculations  are  undertaken. 

2.5  STANDARD  TEST  PROBLEMS. 

As  a  first  standard  test  problem,  the  anomolous  wall  heating  problem 
discussed  above  was  solved  as  the  reflection  of  a  planar  shock  against  a 
vertical  wall  In  2-0.  It  was  found  using  slip  boundary  conditions  that  any 
cut  of  the  2-D  shock  profile  at  constant  -y  replicated  the  1-D  results.  It 
was,  In  fact,  found  that  the  2-0  MFE  solutions  of  this  problem  required  fewer 
time  steps  for  running  to  completion  than  the  1-0  code  required  for  the  cor¬ 
responding  problem  to  be  executed  in  1-0.  The  total  computational  cost  in 
2-0  nevertheless  exceeded  the  cost  for  computing  this  problem  in  1-0,  because 
the  respective  costs  per  time  step  are  much  greater  in  2-D  than  in  1-D.  But 
these  results  suggest  that  the  additional  degrees  of  topological  freedom  for 
MFE  node  movements  in  2-0,  vis  a  vis  1-0,  are  conducive  to  intrinsic  computa¬ 
tional  economies  In  integration  time  steps. 

Physically,  this  benchmark  problem  Is  significant  because  It  contains 
only  fluid  compression  and  expansion  effects;  fluid  shearing  and  vorticlty 
generation  Is  non-existent.  This  test  problem  therefore  acts  as  an  effective 
screening  example  for  both  coding  and  physics  validation  purposes.  It  would 
serve  similarly  as  a  useful  example  for  measuring  numerical  dissipation 

effects  In  all  existing  alrblast  codes. 

The  isodensity  contour  (Isopycnlc)  data  from  shock  tube  measurements  by 
Glass  and  co-workers1  are  perhaps  the  most  highly  pertinent  data  in  existence 
for  the  fundamental  study  of  those  detailed  shock-boundary/ shear  layer  inter¬ 
actions  which  may  ultimately  govern  the  transient  evolution  of  irregular 
shock  reflection  processes.  This  data  resolves  the  density  profiles  in  the 
boundary  and  internal  shear  layers.  The  case  of  regular  shock  reflection  is 
shown  schematically  in  Figure  12. 

The  second  standard  test  considered  in  this  work  is  the  case  of  regular 

shock  reflection  in  argon  at  Mach  nunber  M  =  2.05,  p0  =  150  torr,  p0  = 

3.23  x  IQ"*  g/cm-1,  TQ  =  297. 6°K,  and  wedge  angle  0W  =  60°.  This 


Incident  Shock  ( I ) 


•Reflected  Shock  (R) 


Figure  12.  Schematic  representation  of  regular  reflection  in 
Shock-on-wedge  experiments. 


shock  has  been  measured  experimentally  by  Deschambault  and  Glass.1  Initial 
conditions  are  calculated  from  the  Ranklne-Hugonlot  relationships  In  the 
following  form: 


U  -  v+ 

\  "  ~  C"  “  Ms  for  v+  3  0 


C+  -  (Yp  +  /P+)ls  , 

2ymJ  -  (y  -  1) 
P-  "  P*  *  — (y  ♦  1) - 


rfhere  *  5/3  for  argon. 


(u  -  v  ) 


or  v  * 


P+  -  P. 


P+(U  -  V+)  +  ■  v+)  » 


P-  -  p+ 

P  ♦  U 


(24a) 


(24b) 


p, 


(25) 


*  P+  u  y  • 

For  monatomic  gases,  the  internal  energy  per  unit  volume  Is 
3 

pe  *  215  ‘  (26) 

The  Initial  conditions  are 

P+  -  3.23  x  10-4  g/cm3 
m+  =  0 

E+  *  [P+  e+  +  m£/2P+]  =  3.0  x  10®  ergs/crP  and 
p_  *  7.87  x  10-4  g/cn|3 
m_  *  30.559  gm  -  cm/ sec 
E-  *  [p_  e_  +  m?/2P_]  *  2.132  x  10®  ergs/ cm** 

2.6  BOUNDARY  CONDITIONS. 

A  significant  point  of  fundamental  physics  became  manifest  In  the  study 
of  boundary  conditions  for  shock-on-wedge  calculations.  For  a  variety  of 
reasons,  conventional  alrblast  codes  use  si  Ip- type  boundary  conditions.  In 
actuality,  the  velocity  of  air  molecules  Is  zero  at  the  types  of  surfaces 
under  present  consideration.  For  the  sake  of  common  comparisons,  and  to  test 
the  sensitivity  of  MFE  shock-on-wedge  computations  to  alternative  boundary 
conditions,  we  first  attempted  to  solve  the  regular  shock-on-wedge  test 
problem  with  si  Ip- type  boundary  conditions,  in  which  the  normal  component  of 
velocity  along  the  bottom  horizontal  surface  and  inclined  wedge  are  zero  and 
the  tangential  velocity  components  are  treated  by  zero-Neunann  (or  symmetry) 
conditions  (see  Figure  11).  These  symmetry  boundary  conditions  are  applied 
similarly  to  the  MFE  grid  nodes  so  that  the  grid  nodes  would  slide  freely 
along  both  the  horizontal  surface  and  the  inclined  wedge  surface.  But  the 
MFE  method  has  always  been  found  to  be  extremely  sensitive  In  its  accuracy 
and  consistency  requirements  to  the  presence  of  any  physical  or  mathemati¬ 
cally  Ill-posed  condition.  Here,  the  surface  normal  is  not  defined  uniquely 
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at  the  front  of  the  wedge.  Consequently,  the  basic  consistency  requirements 
in  the  MFE  method  impeded  the  numerical  integration  process  because  the 
symmetry  boundary  condition  at  the  foot  of  the  wedge  is  111 -posed,  and  thus 
could  not  be  properly  resolved  as  a  truly  rigorous  PDE  solution.  The  MFE 
method  also  tends  to  be  less  tolerant  than  most  PDE  methods  of  computational 
swindles  which  are  frequently  used  to  accomplish  slip-type  boundary  condi¬ 
tions.  It  thus  became  quite  apparent  that  It  Is  possible  to  enforce  si  ip- type 
boundary  conditions  only  by  accepting  erroneous  numerical  solutions  In  the 
boundary  layer  region  near  the  front  corner  of  the  wedge.  It  was  further 
apparent  that  these  errors  could  propagate  to  regions  well  away  from  the 
local  source  of  difficulty.  In  view  of  these  results,  we  next  attempted  to 
take  a  giant  step  forward  in  terms  of  both  the  physics  and  computations  of 
the  physically  required  non-slip  boundary  conditions.  This  turned  out  to  be 
a  giant  step  because  this  shock-on-wedge  problem  has  a  turbulent  boundary 
layer  in  the  duct  represented  by  our  problem  domain.  Nevertheless,  some  key 
points  became  manifest  in  this  task  which  will  be  described  immediately 
below. 

The  initial  shock  conditions  are  shown  in  Figure  13.  This  problem  is 
solved  in  a  1  cm  duct.  For  viscosity  y  =  5  x  10"*,  the  Reynolds  number 
is  approximately  5.7  x  104.  Such  a  flow  is  clearly  turbulent.  The  MFE  solu¬ 
tions  of  the  laminar  Navier-Stokes  equations  were  nevertheless  attempted  as  a 
nunerical  experiment.  These  results  appear  in  Figures  14  to  23  with  self- 
explanatory  captions.  At  t  =  0.30  the  shock  is  just  starting  to  encounter 
the  wedge  front;  and  the  Isodensity  contours  early  In  the  reflection  process 
at  t  *  0.35  have  a  regular  profile.  At  t  =  0.44  the  formation  and  subsequent 
shedding  of  eddies  are  becoming  apparent.  These  are  formed  and  evolved 
according  to  the  physically  Inadequate  laminar  dissipation  processes  in  the 
present  Navier-Stokes  solutions  and  not  by  numerical  dissipation  or  by  an 
appropriate  turbulence  model.  Analysis  of  vorticity  generation  rates 
Indicates  clearly  that  a  turbulent  model  of  the  physical  dissipation  in 
the  boundary  layer  is  required  in  order  to  model  this  system  according  to 
physical  reality. 

Mostly,  this  example  shows  an  apparent  capacity  of  the  MFE  method  to 
attempt  to  resolve  extremely  micro-scale  dissipation  processes  simul taneously 
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Figure  13.  Initial  conditions  for  regular  reflection  of  planar 
shock  In  the  experiments  of  Ueschambault  and  Glass. 
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Figure  30.  Projection  of  MFE  Grid  Mesh  on  the  x-y  Plane  for  an 
Experimental  Plane  Shock  of  Deschambault  and  Glass 
Reflecting  Against  a  Vertical  Wall. 

(M$  =  2.05,  8  =  90°,  p+  =  150  torr,  p+  =  3.23  x 
10  ^g/cm^  in  Argon).  A  3  x  31  MFE  Grid  is  Used. 


with  the  macroscopic  flow  features.  In  this  experimental  example,  it  was 
undoubtedly  the  physics  model  which  requires  improvement  more  than  the 
nunerics.  Because  local  turbulent  dissipation  rates  exceed  laminar  rates  by 
large  factors,  it  is  likely  that  this  experiental  example  imposes  more  severe 
demands  upon  the  POE  method  than  will  be  encountered  in  more  realistic  physi¬ 
cal  models  of  real  shock  environments  {which  have  much  larger-than-laminar 
dissipation  rates). 

In  order  to  verify  that  the  complex  behavior  in  this  example  was  truly 
associated  with  vorticlty  generation  in  the  boundary  layer,  this  same  inci¬ 
dent  shock  was  propagated  into,  and  reflected  from,  a  vertical  wall.  Fluid 
shears  do  not  occur  in  this  example,  and  the  MFE  solutions  appear  in  Figures 
24  to  30.  The  shock  profiles  remain  perfectly  planar  (to  many  significant 
figures)  because  only  compressive  and  expansive  forces  act  in  this  reflection 
process.  Neither  numerical  dissipation  nor  anomalous  vorticity  effects  were 
perceptible  in  these  MFE  solutions.  It  was  also  validated  that  the  Rankine- 
Hugonlot  and  Taylor  solutions  of  the  steady  reflected  shock  profiles  were 
obtained  in  these  solutions. 


SECTION  3 
CONCLUSIONS 


The  results  of  work  during  this  phase  of  investigation  continue  to  show 
promise  for  resolving  shock-boundary  layer  interactions  in  2-D  according  to 
PDE's  which  include  the  appropriate  physical  dissipation  effects  of  airblast 
environments.  The  attempted  resolution  of  a  turbulent  shock-boundary  layer 
Interaction  problem  using  only  laminar  dissipation  rates  provided  an  unanti¬ 
cipated  preview  of  physical  modeling  needs  which  may  be  required  in  high- 
resolution  calculations  when  the  nunerical  dissipation  effects  are  success¬ 
fully  reduced  to  insignificant  levels.  In  all  cases,  it  presently  appears 
in  our  MFE  results  that  the  highly  local  physical  dissipation  processes  in 
regions  of  large  fluid  gradients  can  be  sensitive  determinates  of  macroscopic 
flow  properties  in  shock-boundary  layer  interactions.  The  evolution  of  the 
physical  dissipation  effects  over  highly  disparate  scales  in  all  of  the  MFE 
solutions  computed  to  date  exhibit  remarkable  levels  of  agreement  with  those 
corresponding  descriptions  which  have  appeared  in  the  classic  shock  litera¬ 
ture  for  real,  viscous  gases. 


SECTION  4 
RECOMMENDATIONS 


This  project  has  so  far  focussed  on  the  feasibility  of  the  MFE  method 
contributing  to  the  technology  base  for  airblast  studies.  Accordingly,  the 
work  to  date  has  revealed  numerous  research  needs  which  should  now  be 
addressed  in  order  to  push  onward  toward  real-world  applications.  These 
incl ude: 

1.  Grid  restructuring.  The  grid  connections  used  in  the  existing  MFE 
airblast  code  in  2-D  contain  certain  biasing  effects  which  would  be 
undesirable  in  complex  flow  fields.  Alternative  grid  connection 
schemes  are  available  and  should  be  developed  for  DNA  applications. 

2.  Implement  ADI  methods  of  matrix  solution  for  more  efficient  solution 
of  problems  on  large  grid  meshes. 

3.  Extended  benchmarking  of  MFE  code  against  classic  laminar  flows. 
In  the  interest  of  feasibility  projections,  the  testing  of  the 
current  MFE  code  version  has  been  applied  to  preview  nunerical 
requirements  in  some  shock-boundary  layer  regimes  which  are  well 
beyond  the  laminar  Navier-Stokes  model  which  is  contained  in  this 
code  version.  It  is  now  timely  to  drop  back  from  this  advanced 
stage  of  experimentation  and  systematically  perform  those  numerous 
validations  against  some  of  the  basic  laminar  boundary  layer  bench¬ 
mark  problems.  Particular  emphasis  should  be  directed  toward  view¬ 
ing  the  onset  and  development  of  both  laminar  and  possible  early 
turbulent  flow  transitions. 
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